rm(list=ls())
library('dplyr')
library('tidyr')
library('VennDiagram')
library('grid')
library('futile.logger')
library("ggplot2")
library("gplots")
library('plotly')
library('tools')
library('cowplot')
library('UpSetR')
Import des 3 vcf obtenus après lancement du pipeline (vcf_12x, vcf_24x, vcf40x) et import du vcf de référence (vcf_ref).
vcf_12x = read.table("/home/elodie/Documents/Analysis/gold_12x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_24x = read.table("/home/elodie/Documents/Analysis/gold_24x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_40x = read.table("/home/elodie/Documents/Analysis/gold_40x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_ref = read.table("/home/elodie/Documents/Elodie/Datas_ismael_ref/gold_on_data_elodie_ismael_decomposed_normalize_uniq.vcf")
plot1 = ggplot(nombre_variants, aes(fill=type, y=vcf, x=nombre)) +
geom_bar(position = "stack", stat = "identity") +
scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
labs(fill = "Type de variants", x = "Nombre de variants", y = "VCF") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5)) +
labs(title = "Nombre et type de variant par VCF")
ggplotly(plot1)
La distribution des types de variants semble équivalente pour chaque VCF. Cependant dans les 3 VCF générés par le pipeline on trouve des variants “Undetermined” mais pas dans le VCF de référence.
Que sont les Undetermined non retrouvés dans le VCF de référence ?
undetermined = vcf_all %>% filter(type == "undetermined")
print(undetermined)
## variants chromosome position
## 1 chr4_66388036_A_* chr4 66388036
## 2 chr4_66440407_CATATATATACATATATATACAT_* chr4 66440407
## 3 chr4_66214417_A_* chr4 66214417
## 4 chr4_66440407_CATATATATACATATATATACAT_* chr4 66440407
## 5 chr4_66473213_A_* chr4 66473213
## ref alt format Id gene
## 1 A * 1/.:0,5:6:41:205,41,45 423 EPHA5
## 2 CATATATATACATATATATACAT * 1/.:0,16:22:99:867,243,285 535 EPHA5
## 3 A * 1/.:0,11:15:99:717,152,110 110 EPHA5
## 4 CATATATATACATATATATACAT * 1/.:0,22:31:99:1233,367,435 547 EPHA5
## 5 A * 1/.:0,4:15:99:1163,469,449 651 EPHA5
## type vcf
## 1 undetermined vcf_12x
## 2 undetermined vcf_24x
## 3 undetermined vcf_40x
## 4 undetermined vcf_40x
## 5 undetermined vcf_40x
Variants avec * en alt lié au vt décompose ?
pos66388036 = vcf_all %>% filter(variants == "chr4_66388034_GTA_G")
print(pos66388036)
## variants chromosome position ref alt format
## 1 chr4_66388034_GTA_G chr4 66388034 GTA G 0/1:2,5:7:46:162,0,46
## 2 chr4_66388034_GTA_G chr4 66388034 GTA G 0/1:7,8:15:99:250,0,267
## 3 chr4_66388034_GTA_G chr4 66388034 GTA G 0/1:13,15:28:99:472,0,494
## 4 chr4_66388034_GTA_G chr4 66388034 GTA G 0/1:.:362:146,111:176,136:155
## Id gene type vcf
## 1 422 EPHA5 deletion vcf_12x
## 2 441 EPHA5 deletion vcf_24x
## 3 453 EPHA5 deletion vcf_40x
## 4 413 EPHA5 deletion vcf_ref
pos66388036bis = vcf_all %>% filter(variants == "chr4_66388036_A_G")
print(pos66388036bis)
## variants chromosome position ref alt format Id
## 1 chr4_66388036_A_G chr4 66388036 A G ./1:0,1:6:41:205,176,204 424
## gene type vcf
## 1 EPHA5 SNV vcf_12x
A cette position délétion de 2 bases (66388035+66388036) bien détectée dans les 4 VCF Substitution de A>G en 66388036 dans le VCF_12x (comme peu de reads, variant caller se trompe? Faux positif?)
pos66440407 = vcf_all %>% filter(position == 66440407)
print(pos66440407)
## variants chromosome position
## 1 chr4_66440407_CATATATATACATATATATACAT_* chr4 66440407
## 2 chr4_66440407_CATATATATACATATATATACAT_C chr4 66440407
## 3 chr4_66440407_CATATATATACATATATATACAT_* chr4 66440407
## 4 chr4_66440407_CATATATATACATATATATACAT_C chr4 66440407
## ref alt format Id gene
## 1 CATATATATACATATATATACAT * 1/.:0,16:22:99:867,243,285 535 EPHA5
## 2 CATATATATACATATATATACAT C ./1:0,6:22:99:867,631,654 536 EPHA5
## 3 CATATATATACATATATATACAT * 1/.:0,22:31:99:1233,367,435 547 EPHA5
## 4 CATATATATACATATATATACAT C ./1:0,9:31:99:1233,875,889 548 EPHA5
## type vcf
## 1 undetermined vcf_24x
## 2 deletion vcf_24x
## 3 undetermined vcf_40x
## 4 deletion vcf_40x
Zone avec des reads ayant des dététions différentes qui se chevauchent dans les 4 BAM. Le VCF de référence n’appelle pas de variants dans cette zone. Zone difficile pour le variant caller, faire appel à un variant caller spécialiste des délétions?
pos66214417 = vcf_all %>% filter (position == 66214417)
print(pos66214417)
## variants chromosome position ref alt format Id
## 1 chr4_66214417_A_* chr4 66214417 A * 1/.:0,11:15:99:717,152,110 110
## 2 chr4_66214417_A_T chr4 66214417 A T ./1:0,4:15:99:717,468,450 111
## gene type vcf
## 1 EPHA5 undetermined vcf_40x
## 2 EPHA5 SNV vcf_40x
pos66214413 = vcf_all %>% filter (position == 66214413)
print(pos66214413)
## variants chromosome position ref alt format
## 1 chr4_66214413_TATAA_T chr4 66214413 TATAA T 0/1:1,9:10:15:375,0,15
## 2 chr4_66214413_TATAA_T chr4 66214413 TATAA T 0/1:4,11:15:99:450,0,110
## Id gene type vcf
## 1 108 EPHA5 deletion vcf_24x
## 2 109 EPHA5 deletion vcf_40x
Délétion de 4 bases de 66214413 à 66214417 vu dans les 4 BAM mais uniquement appelée dans VCF 24x et 40x. Non appelé dans le VCF de référence. Subtitution en 66214417 aussi un A>T vu dans les 4 BAM mais uniquement appelée dans le VCF 40x. Non appelé dans le VCF référence.
pos66473213 = vcf_all %>% filter (position == 66473213)
print(pos66473213)
## variants chromosome position ref alt format Id
## 1 chr4_66473213_A_C chr4 66473213 A C 1/1:0,18:18:54:720,54,0 635
## 2 chr4_66473213_A_* chr4 66473213 A * 1/.:0,4:15:99:1163,469,449 651
## 3 chr4_66473213_A_C chr4 66473213 A C ./1:0,11:15:99:1163,206,135 652
## gene type vcf
## 1 EPHA5 SNV vcf_24x
## 2 EPHA5 undetermined vcf_40x
## 3 EPHA5 SNV vcf_40x
Dans le VCF de référence en position 66473209 à 66473213 : délétion de 4 bases appelée mais peu de reads avec la délétion et beaucoup plus de reads avec substitutions de A>C en 66473213. Substitutions appelée uniquement dans VCF 24x et 40x. Dans les zones où il y a des délétions et substitutions, le variant caller rencontre des difficultés?
plot2 = ggplot(nombre_variants_par_gene, aes(fill=gene, y=vcf, x=nombre)) +
geom_bar(position="stack", stat="identity") +
scale_fill_manual(values = c("#fdc086", "#91bfdb", "#b2df8a")) +
labs(fill = "Nom des gènes", x = "Nombre de variants", y = "VCF") +
theme(panel.grid = element_blank(),plot.title = element_text(hjust = 0.5)) +
labs(title = "Nombre de variants par gène et par VCF")
ggplotly(plot2)
Le nombre de variants est identique dans chaque VCF pour les gènes UTS2 et LPL mais le nombre varie pour le plus grand gène UPHA5.
nombre_variants_par_gene$vcf_nom_genes = paste(nombre_variants_par_gene$vcf, nombre_variants_par_gene$gene, sep="_")
plot3 = ggplot(nombre_variants_par_gene, aes(fill=type, y=vcf_nom_genes, x=nombre)) +
geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
labs(title = "Nombre de variants par gène, par VCF et par type de variant") +
labs(fill = "Type de variants", x = "Nombre de variants", y = "VCF") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
ggplotly(plot3)
La répartition des types de variants par gène et par VCF est concordante. Aucune insertion dans UTS2 dans les 4 VCF.
plot4 = ggplot(nombre_type_snv[!is.na(nombre_type_snv$type), ], aes(fill = type, y = vcf, x = nombre))+
geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
scale_fill_manual(values = c("#fc9272", "#9ebcda")) +
labs(title = "Type de SNV par VCF", x = "Type de substitutions", y = "Nombre de substitutions") +
labs(fill = "Type de SNV") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
plot4
En génétique humaine, il y a un environ 2 fois plus de transitions que de transversions. Sur le plot, on voit environ 2/3 de transitions et 1/3 de transversions ce qui est cohérent. Valeurs cohérentes dans tous les VCF.
plot5 = ggplot(nombre_type_snv[!is.na(nombre_type_snv$type_snv), ], aes(fill = type_snv, y = vcf_gene, x = nombre)) +
geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
labs(title = "Répartition des types de SNV par VCF et par gène", x = "Nombre", y = "VCF et gène") +
scale_fill_manual(values = c("#a6cee3", "#1f78b4","#b2df8a", "#33a02c", "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a", "#ffff99", "#b15928")) +
labs(fill = "Type de SNV", x = "Nombre", y = "VCF") +
theme(plot.title = element_text(hjust = 0.5, size=14, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
theme(panel.grid = element_blank())
ggplotly(plot5)
Une répartition des types de SNV cohérente entre VCF pour chaque gène.
plot6 = ggplot(vcf_all %>% filter (type == "insertion" | type == "deletion"), aes(x = taille_delins, y = vcf)) +
geom_boxplot(aes(fill=type)) +
theme_classic() +
labs(x = "Taille des delins", y = "VCF", title = "Boxplots de la taille des insertions et délétions par VCF") +
theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
labs(fill = "Type de delins") +
scale_fill_manual(values = c("deletion" = "#bebada", "insertion" ="#ffffb3"))
plot6
La distribution des tailles des délétions est équivalente entre les 4 VCF. Les VCF 12x, 24x et 40x contiennent plus de valeurs aberrantes pour les insertions et délétions que le VCF de référence. La taille maximale des insertions du VCF de référence est plus petite que pour les 3 VCF obtenus avec le pipeline. L’interquartile est plus petit également pour la taille des insertions du VCF de référence ce qui signifie que les valeurs des insertions sont plus regroupées autour de la médiane. Il y a moins de disparité des valeurs pour le VCF de référence.
plot7 = ggplot(vcf_all %>% filter (type == "insertion"), aes(x = taille_delins, y = vcf)) +
geom_boxplot(aes(fill=gene)) +
theme_classic() +
labs(x = "Taille des insertions", y = "VCF", title = "Boxplots de la taille des insertions par VCF et par gène") +
theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
labs(fill = "Noms des gènes") +
scale_fill_manual(values = c("EPHA5" = "#fdc086", "LPL" ="#91bfdb", "UTS2" = "#b2df8a"))
plot8 = ggplot(vcf_all %>% filter (type == "deletion"), aes(x = taille_delins, y = vcf)) +
geom_boxplot(aes(fill=gene)) +
theme_classic() +
labs(x = "Taille des délétions", y = "VCF", title = "Boxplots de la taille des délétions par VCF et par gène") +
theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
labs(fill = "Noms des gènes") +
scale_fill_manual(values = c("EPHA5" = "#fdc086", "LPL" ="#91bfdb", "UTS2" = "#b2df8a"))
plot_grid(plot7, plot8)
plot9 = ggplot(vcf_12x_24x_40x, aes(x = as.numeric(vcf_12x_24x_40x$DP), y = vcf)) + geom_boxplot(aes(fill=gene)) + theme_classic() +
theme_classic() +
labs(x = "Profondeur de lecture des variants", y = "VCF", title = "Boxplots de la profondeur de lecture des variants par gène") +
theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
labs(fill = "Noms des gènes") +
scale_fill_manual(values = c("UTS2" = "#b2df8a", "LPL" ="#91bfdb", "EPHA5" = "#fdc086"))
plot9
## Warning: Use of `vcf_12x_24x_40x$DP` is discouraged.
## ℹ Use `DP` instead.
plot10 = ggplot(dp_12x, aes(x = DP, y = nombre)) +
geom_histogram(stat = "identity", fill = "#a6cee3") +
geom_vline(xintercept = mean(dp_12x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 12x") +
theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5, size = 14, face = "bold")) +
theme(panel.grid = element_blank())
## Warning in geom_histogram(stat = "identity", fill = "#a6cee3"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot11 = ggplot(dp_24x, aes(x = DP, y = nombre)) +
geom_histogram(stat = "identity", fill = "#1d91c0") +
geom_vline(xintercept = mean(dp_24x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 24x") +
theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5, size = 14, face = "bold")) +
theme(panel.grid = element_blank())
## Warning in geom_histogram(stat = "identity", fill = "#1d91c0"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`
plot12 = ggplot(dp_40x, aes(x = DP, y = nombre)) +
geom_histogram(stat = "identity", fill = "#253494") +
geom_vline(xintercept = mean(dp_40x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 40x") +
theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5) ) +
theme(panel.grid = element_blank())
## Warning in geom_histogram(stat = "identity", fill = "#253494"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`
plot_grid(plot10, plot11, plot12, ncol = 3, nrow = 1)
On ne fera pas d’exploration de la VAF ni de la profondeur DP pour le VCF de référence car il n’y a pas d’explication aux valeurs de DP, AD, ADALL dans le VCF. Les valeurs ne sont pas cohérentes avec ce qui est observé dans le bam de référence associé au VCF dans IGV.
cat("VAF minimale dans VCF 12x :", min(vaf_12x$vaf), "%", "\n")
## VAF minimale dans VCF 12x : 16.7 %
cat("VAF minimale dans VCF 24x :", min(vaf_24x$vaf), "%", "\n")
## VAF minimale dans VCF 24x : 14.3 %
cat("VAF minimale dans VCF 40x :", min(vaf_40x$vaf), "%", "\n")
## VAF minimale dans VCF 40x : 16.3 %
plot13 = ggplot(vcf_vaf, aes(x = nombre_variants, y = cat_vaf , fill = vcf)) +
geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
labs(title = "Répartition du nombre de variants par catégorie de VAF", x = "Nombre de variants", y = "Catégorie de VAF") +
scale_fill_manual(values = c("#a6cee3", "#1d91c0", "#253494")) +
labs(fill = "VCF") +
theme(plot.title = element_text(hjust = 0.5, size=14, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
theme(panel.grid = element_blank()) +
theme_classic()
ggplotly(plot13)
variants_list = list(vcf_12x = unique(vcf_12x$variants), vcf_24x = unique(vcf_24x$variants), vcf_40x = unique(vcf_40x$variants), vcf_ref = unique(vcf_ref$variants))
plot14 = upset(fromList(variants_list), order.by = "freq", text.scale = 2, point.size = 5, sets.bar.color = "#999999", main.bar.color = "#fc8d62")
plot14
plot15 = display_venn(variants_list, category.names = c("vcf_12x" , "vcf_24x" , "vcf_40x", "vcf_ref"),
fill = c("#a6cee3", "#1d91c0", "#253494", "#E69F00"), cat.cex = 1.2,
cat.fontface = "bold",
cat.default.pos = "outer",
cat.dist = c(0.1, 0.1, 0.1, 0.1))
plot16 = ggplot(cat_vcf_all_nombre, aes(fill=type, x=nombre, y=categorie)) +
geom_bar(position = "stack", stat="identity") +
scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
labs(title = "Nombre de variants par catégories et par type de variants") +
labs(fill = "Type de variants", x = "Nombre de variants", y = "Catégories de VCF") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
plot17 = ggplot(cat_vcf_all_genes, aes(fill=gene, x=nombre, y=categorie)) +
geom_bar(position = "stack", stat="identity") +
scale_fill_manual(values = c("UTS2" = "#b2df8a", "LPL" ="#91bfdb", "EPHA5" = "#fdc086")) +
labs(title = "Nombre de variants par catégories et par gènes") +
labs(fill = "Type de variants", x = "Nombre de variants", y = "Noms des gènes") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
plot_grid(plot16, plot17, ncol = 1, nrow = 2)
variants_list_12x_24x_40x = list(vcf_12x = unique(vcf_12x$variants), vcf_24x = unique(vcf_24x$variants), vcf_40x = unique(vcf_40x$variants))
plot18 = display_venn(variants_list_12x_24x_40x, category.names = c("vcf_12x" , "vcf_40x", "vcf_24x"),
fill = c("#a6cee3", "#1d91c0", "#253494"), cat.cex = 1.2,
cat.fontface = "bold",
cat.default.pos = "outer",
cat.dist = c(0.1, 0.1, 0.1))
recall (sensibilité), précision (spécificité) et score F1. Ce sont des métriques qui permettent de voir si un modèle est performant. Le score F1 combien recall et précision Plus le recall est élevé, plus le modèle repère des positifs. Plus la précision est élevée, moins le modèle se trompe sur les positifs (moins il y a de faux positifs). Plus le score F1 est élevé, plus le modèle est performant.
happy_12x_summury = read.table("/home/elodie/Documents/Analysis/Happy_12x/Happy_12x/ref-12x.summary.csv", header=TRUE, sep=',')
happy_24x_summury = read.table("/home/elodie/Documents/Analysis/Happy_24x/Happy_24x/ref-24x.summary.csv", header=TRUE, sep=',')
happy_40x_summury = read.table("/home/elodie/Documents/Analysis/Happy_40x/Happy_40x/ref-40x.summary.csv", header=TRUE, sep=',')
plot_roc_snp = plot_data(roc_snp, is.PR=TRUE)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot_roc_indel = plot_data(roc_indel, is.PR=TRUE)
plot_grid(plot_roc_snp, plot_roc_indel, nrow=2, ncol=1)